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In searches for gravitational waves emitted by known isolated pulsars in data collected by a 
detector one can assume that the frequency of the wave, its spindown parameters, and the position 
of the source in the sky are known, so the almost monochromatic gravitational-wave signal we are 
looking for depends on at most four parameters: overall amplitude, initial phase, polarization angle, 
and inclination angle of the pulsar's rotation axis with respect to the line of sight. We derive two 
statistics by means of which one can test whether data contains such gravitational- wave signal: the 
C/-statistic for signals which depend on only two unknown parameters (overall amplitude and initial 
\ phase), and the J^-statistic for signals depending on all four parameters. We study, by means of 

^Sj . the Fisher matrix, the theoretical accuracy of the maximum-likelihood estimators of the signal's 

parameters and we present the results of the Monte Carlo simulations we performed to test the 
accuracy of these estimators. 
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I. INTRODUCTION 



qh' We study the detection of almost monochromatic gravitational waves emitted by known single pulsars in data 
I . collected by a detector. Several such searches were already performed with data collected by the LIGO and GEO600 
^jPj' detectors [iHsj. We thus assume that the frequency of the wave (together with its time derivatives, i.e. the spindown 
\ parameters) and the position of the source in the sky are known. The gravitational- wave signal we are looking for 
depends on at most four (often called amplitude) parameters: overall amplitude, initial phase, polarization angle, and 
inclination angle (of the pulsar's rotation axis with respect to the line of sight). 
^ , In Sec. 2 we introduce three statistics by means of which one can test whether data contains a gravitational- 
' wave signal: the "H-statistic for completely known signals, the t?-statistic for signals which depend on only two 
unknown parameters (overall amplitude and initial phase), and the J^-statistic suitable for signals depending on all 
four amplitude parameters. Both statistics Q and J- are derived from the maximum likelihood (ML) principle, and 
, the statistic Q is independently obtained using Bayesian approach and the composite hypothesis testing. In Sec. 3 we 
' study, by means of the Fisher matrix, the theoretical accuracy of the ML estimators of the signal's parameters and in 
\ Sec. 4 we present the results of the Monte Carlo simulations we performed to test the accuracy of the ML estimators. 



II. USING THE F AND g STATISTICS TO PERFORM TARGETED SEARCHES FOR 
GRAVITATIONAL WAVES FROM PULSARS 

In the case when the signal s{t) we are looking for is completely known, the test that maximizes probability of 
detection subject to a certain false alarm probability is the likelihood-ratio test, i.e. we accept the hypothesis that 
the signal is present in detector's data x if 

A(x):=^>Ao, (1) 

where the likelihood function A(a;) is the ratio of probability densities pi{x) and Pq{x) of the data x when the signal 
is respectively present or absent. The parameter Aq is a threshold calculated from a chosen false alarm probability. 
Assuming stationary and additive Gaussian noise with one-sided spectral density constant (and equal to Sq) over the 
bandwidth of the signal, the log likelihood function is approximately given by [y] 

lnA[a:(t)] - 2^ [{x{t)s{t)) - \{s{tf)^ , (2) 

where To is the observation time and the time- averaging operator (•) is defined as 

(5):=^ r°9it)dt. (3) 
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Equation ([2]) implies that the likehhood-ratio test ([T]) can be replaced by the test 

n[xit)] {x{t)s{t)) > Ho, (4) 

where the optimal statistic H in this case is the matched filter and Ho is the threshold for detection. 

Suppose now that the signal s{t; 6) depends on a set of unknown parameters 6, then a suitable test can be obtained 
using a Bayesian approach and composite hypothesis testing. The composite hypothesis in this case is the hypothesis 
that when a signal is present it can assume any values of the parameters. Assuming that the cost functions are 
independent of the values of the parameters, we obtain the following Bayesian decision rule to choose the hypothesis 
that the signal is present (see e.g. 0, Chapter 5.9): 

— ^ / Pi(a;;0)7r(0)d0>7o, (5) 
Poix) Je 

where Q is the parameter space on which 9 is defined and Tr{6) is the joint a priori distribution of 9. The expression 
on the left hand side of Eq. ^ is know as the Bayes factor and it is the ratio between the posterior probability 
distribution on the signal parameters marginalized over the parameters themselves (this is the signal model Bayesian 
evidence) and the noise model which has no defining parameters (this is the noise model Bayesian evidence). 

As a template for the response of an interferometric detector to the gravitational- wave signal from a rotating neutron 
star we use the model derived in @. This template depends on the set of following parameters: 9 = {ho, 4>o,i^, l, f , S, a), 
where ho is the dimensionless amplitude, 0o is an initial phase, ip is the polarization angle, t is the inclination angle, 
angles d (declination) and a (right ascension) are equatorial coordinates determining the position of the source in 
the sky, and the 'frequency vector' f := (/o, /i, /2, • . • ) collects the frequency fo and the spindown parameters of 
the signal. In the case of pulsars known from radio observations we in general know the subset ^ — (f , S, a) of the 
parameters 9. 

Sometimes, like in the case of the Vela pulsar, we also know from X-ray observations the values of the angles ip and 
I (see 0,0 for observational results). We then have only two unknown parameters: ho and (f>o- The response s{t) of 
the detector to the gravitational wave we can write in this case in the following form ^] : 

s{t) — ho cos (f)o hc{t) -\- ho sin^o hgit), (6) 

where he and hs are known functions of time, 

hJt) := A+{cofi2%phi{t) + sin2V'/i2(0) - Ax ( sin 2^- /i3(i) - cos 2iIj hAt)) , 

(7) 

h,{t) -Ax(sin2'0/ii(i) - cos 2?/; /i2(t)) - A+( cos 2V' /i3(t) -t- sin2V' /i4(0) • 
Here the constants A+ and Ax are 

A+ := i(l + cos^ t), Ax := cost, (8) 

and the four functions of time /i^ (fc = 1, . . . , 4) depend only on parameters ^ and are defined as follows 
hi{t; $,) := a{t; 5, a) cos (j){t; f , S, a), h2{t; ^) := b{t; S, a) cos <j){t; f , S, a), 

(9) 

hsit; ^) := a{t; S, a) sin (f>{t; f , 6, a), hi{t; ^) := b{t; 5, a) sin (j}{t; f , 5, a), 

where a, b are the amplitude modulation functions and (j) is the phase modulation function. Their explicit forms are 
given in Q. 

Let us calculate the likelihood function for the signal ([6]). Observing that the amplitude modulation functions a 
and b vary much more slowly than the phase 4> of the signal and assuming that the observation time is much longer 
than the period of the signal we approximately have 

(/ll hs) ^ (hi hi) ^ {h2 hs) 9^ {h2 hi) 0, 

(/ll h{) = {hs hs) = \A, {h2 h2) = (/14 hi) - IB, {hi h2) = {hs hi) = \G, 

where we have introduced the time averages 

A:=(a2), B:={b^), C:^{ab). (11) 

As a consequence of the above approximations we have the following approximate expressions for the time averaged 
products of the functions he and hg, 

{hl)^{hl)^N, (W=0, (12) 
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where TV is a constant defined as 

TV := i (^A{Al_ cos^ 2?/' + Al sin^ 2?/') + B{Al sin^ 2V^ + cos^ 2^-) 

+ C(A^- A2)sin4V'). (13) 
With the above approximations the hkehhood function A for the signal ([6|) can be written as 

lnA[x{t)-(j)o,ho] ^ 2^ (^hocosc^o{x{t)h,{t)) + hosinc^o{x{t)h,{t)) - ^hlN^ . (14) 

Let us also note that the optimal signal-to-noise ratio (SNR) p for the signal ^ (see [s] for definition) can be 
approximately computed as 



It is natural to assume that the prior probability density of the phase parameter 0o is uniform over the interval 
[0, 27r) and that it is independent of the distribution of the amplitude parameter Hq, i.e. 

7r(0o) = ^, 0e[O,2^). (16) 

With the above assumptions the integral J^^ pi{x; (/)q, /io)7r(0o)d0o can be explicitly calculated (see [9], Chapter 7.2) 
and we obtain the following decision criterion 



exp ( - lo ( 2hJ^g[xit)] ] > 70, (17) 



JQ / \ V -JO 

where /q is the modified Bessel function of zero order and the statistic G is defined as 

G[xit)] := J^^i^{xit)h,it)f + {x{t)h,it)f)- (18) 

The function on the left-hand side of Eq. (jl7p is a monotonically increasing function of Q and it can be maximized if 
G is maximized independently of the value of Hq. Thus the test 

G[xit)] > Go, (19) 

provides a uniformly most powerful test with respect to the amplitude /iq- 

When we have no a priori information about the parameters a standard method is the maximum likelihood (ML) 
detection which consists of maximizing the likelihood function A[x{t); 9] with respect to the parameters of the signal. 
If the maximum of A exceeds a certain threshold we say that the signal is detected. The values of the parameters 
that maximize A are said to be the ML estimators of the parameters of the signal. For the case of signal ^ it is 
convenient to introduce new parameters 

Ac ■— /iQ cos As :— /iQsin^o- (20) 
Then one can find the ML estimators of the amplitudes Ac and A^ in a closed analytic form, 

Ac- ^ , As- ^ . (21) 

It is easy to find that the estimators Ac and A^ are unbiased and also that they are of minimum variance, i.e. their 
variances attain the lower Cramer- Rao bound determined by the Fisher matrix. The variances of both estimators are 
the same and equal to 1/A'^. Substituting the estimators Ac and As for the parameters Ac and As in the likelihood 
function one obtains a reduced likelihood function. This reduced likelihood function is precisely equal to the G- 
statistic given by Eq. (|18p . i.e. tj[a::(i)] — In A[a;(<); ^c, ^s]- The formula for the CJ-statistic obtained without usage of 
the simplifying assumptions (llOp is given in Appendix A. 

When the all four parameters (/iq, 00, '-) are unknown one can introduce new parameters Ak [k — 1, . . . , 4) that 
are functions of (/iq, ^o, V', '-) such that the response s{t) takes the form 

s{t)=Aihi{t) + A2h2{t) + A3h3{t)+Aihi{t), (22) 
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FIG. 1: Receiver operating characteristic (ROC) for the statistics "H, Q, and T for the optimal signal-to-noise ratio p = 2. 



where the functions hk are given by Eqs. © and the parameters read 

Ai /io+ cos 2'0 cos 0o — /lox sin 2i/; sin 
A2 := /10+ sin 2?/; cos 0o + ^ox cos 2-0 sin 

(23) 

A3 : = — /10+ cos 2t]j sin (po — hox sin 2^ cos 0o , 
A4 := — /io+ sin 2?/' sin (/io + ^ox cos 2?/; cos (/)o; 

here /10+ := '^0^4+ and ft-ox ^o^x [see Eq. ([5])]. The ML estimators of Ak can again be obtained in an explicit 
analytic form and the reduced likelihood function is the J^-statistic given by (see Q for details) 

2T / 

J^[x{t)] := In A[x(i); ii, . . . , i4] = ^ {{xhif + (xhsf) + A {{xh^f + {xh^f) 

- 2C {{xhi){xh2) + {xh3){xh4)y (24) 

where D := AB -C'^. The test 

T[xit)] > To (25) 

is not a uniformly most powerful test with respect to unknown parameters {ho, (/>o, V': '-)• It was recently shown that 
uniform a priori distributions of (/iq, 0Oi V"; cos <.) lead to a statistic that can be more powerful than T [11] . 

In Fig.[l]we have plotted the receiver operating characteristics (ROC) for the three statistics H, G, and T considered 
in the present section. 

III. THE FISHER MATRIX 

Using the Fisher matrix we can assess the accuracy of the parameter estimators. We have two theorems that can 
loosely be stated as follows. 

Theorem 1 (Cramer-Rao bound) The diagonal elements of the inverse of the Fisher matrix are lower bounds on 
the variances of unbiased estimators of the parameters. 

Theorem 2 Asymptotically (i.e. when the SNR tends to infinity) the ML estimators are unbiased and their covariance 
matrix is equal to the inverse of the Fisher matrix. 
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FIG. 2: Dependence of standard deviations (calculated from the Fisher matrix) of the parameters ho, (j>o, ip, and cost on the 
cosine of the inclination angle t. We have taken (po ~ 4.03 and ^ = —0.22 (values of other parameters needed to perform the 
computation of the Fisher matrix are listed in the text of Sec. 3). 



For an almost monochromatic signal s ~ s{t] 6), which depends on the parameters 9 = (9i, . . . , 0„i), the elements 
of the Fisher matrix F can be approximately calculated from the formula 

Te.e, = Tj^Tj^ ), i,j = l,...,m. (26) 



In the case when only the parameters /iq and are unknown (t/-statistic search), the Fisher matrix can be computed 
easily from Eqs. (jH) and (|26|) . It is diagonal and the standard deviations of the parameters defined as the square 
roots of the diagonal elements of the inverse of the Fisher matrix read: 

ho p p 

where p is the optimal SNR [given in Eq. ([T5|) ]. 

When all the four amplitude parameters ho, (f>o, "0, and l are unknown (J-"- statistic search), the Fisher matrix can 
be computed by means of formulas given in Appendix B. In this case it is not diagonal, indicating that the amplitude 
parameters are correlated. The quantities au^/ho, cr^o, t^, (where the standard deviations again are defined as 
square roots of diagonal elements of the inverse of the Fisher matrix) have rather complicated analytical form but 
they possess a number of simple properties. They are inversely proportional to the overall amplitude ho, independent 
on the initial phase 0O: and very weakly dependent on ?/;, however there is a strong dependence on l. 

In Fig. [2] we have shown the dependence of the standard deviations on the cosine of the inclination angle t. The 
time averages from Eqs. (needed to compute the Fisher matrix) were computed here for the location of the Virgo 
detector [IJl and for a randomly chosen position of the source in the sky. We have also taken ho — 6.0948 x 10^^, 
To = 441610 s, and S'o = 2 Hz^i, which corresponds to the SNR p = 28MV2N [see Eq. ([15])]. The same time 
averages and the values of To, ho, Sq were used in the Monte Carlo simulations described in Sec. 4. We see in Fig. 
[2] that the standard deviations become singular when cost = ±1. This singularity originates from the degeneracy of 
the amplitude parameters for cost = ±1. In this case the amplitude parameters from Eqs. (j23p become 



Ai = /ioCOs(2V'±0o), ^2 = /iosin(2V'±0o), A3 = T^2, A4 = ±Ai. (28) 
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cos I = 0.1 cos I = -0.93 




FIG. 3: Mean and normalized standard deviation of the ML estimator of the amplitude ho as a function of the SNR. The top 
two panels are the means of the estimator for the two values of cos t. The continuous line is the true value and the circles 
are results of the simulation for 1000 realizations of the noise. The bottom two panels are the standard deviations. The 
continuous line is obtained form the Fisher matrix whereas the circles are results of the simulation. We have taken (po ~ 4.03 
and Jp = -0.22. 



Thus only two of them arc independent. Therefore the determinant of the 4-dimensional Fisher matrix is equal to 
zero at cosi = ±1 and consequently its inverse does not exist in this case. 

IV. MONTE CARLO SIMULATIONS 

We have performed two Monte Carlo simulations in order to test the performance of the ML estimators. We have 
compared the simulated standard deviations of the estimators with the ones obtained from the Fisher matrix. In 
particular we have investigated the behavior of the ML estimators near the Fisher matrix singularity at cost = ±1. 
In each simulation run we have generated the signal using Eq. (j22p . we have added it to a white Gaussian noise, and 
we have estimated the amplitude parameters using the J^-statistic. Each simulation run was repeated 1000 times for 
different realizations of the noise. 

In the first simulation we have investigated the bias and the standard deviation of the ML estimator of the amplitude 
parameter ho as functions of the SNR for the two cases: cost = 0.1 and cost = —0.93. The results are presented in 
Fig. |3l For the first case the ML estimator is nearly unbiased and its standard deviation is close to the one predicted 
by the Fisher matrix even for low SNRs. In the second case the simulation shows considerable bias of the estimator 
and its standard deviation lower than the one predicted by the Fisher matrix. However, Theorem 2 is satisfied in the 
second case. For cos t close to ±1 we have to go to SNR ^ 1000 in order for the ML estimator to be unbiased and its 
standard deviation close to the one given by the Fisher matrix. 

In the second simulation, illustrated in Fig. 4, we have investigated the bias and the standard deviation of the ML 
estimators of the amplitude parameters Hq and cost as functions of cost for the fixed SNR p = 15.6. We find that 
for I cost| < 0.5 the biases are less than 10% and the Fisher matrix overestimates the standard deviations also by less 
than 10%. We see that over the whole range of cost the standard deviations of the parameters are roughly constant 
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FIG. 4: Means and normalized standard deviations of the ML estimators of ho and cos t as functions of cos t. The top two 
panels are the means of the estimators. The continuous lines are the true values and the circles are results of the simulation for 
1000 realizations of the noise. The bottom two panels are the standard deviations. The continuous lines are obtained form the 
Fisher matrix whereas the circles are results of the simulation. We have assumed (po ~ 4.03, ip = —0.22, and p — 15.6. Plots 
for < cost < +1 (not shown here) are mirror images of the plots for ~1 < cost < 0. 



whereas the biases increases as the | cos l\ increases. At cos til the amplitude Hq is overestimated by ahnost a factor 
of 2. 

One reason why Theorem 1 does not apply here is that it holds for unbiased estimators. Also a more precise 
statement of Theorem 1 (see e.g. Theorem 8 in [lo| ) requires that the Fisher matrix F is positive definite for all values 
of parameters. This last assumption is clearly not satisfied here as det F for cos t = ±1. 



Appendix A: The general form of the (/-statistic 



It is not difficult to obtain the tj-statistic without simplifying assumptions (jlOp . The estimators of the amplitude 
parameters Ac and Ag are then given by 

^ ^ {hl){xhc) - {hchs){xhs) ^ ^ {hl){xh,) - {hch,){xhc) 
{hl){h^ ^ {Khsf ' ' {hl){h^) - {hchsf ' 

and the general form of the CJ-statistic reads 



( {hl){xKf -2{Kh,){xK){xhs) + {hl){xKY 

sA {hi){hi)-{hcKY 



Gi^it)] = — ( ^ 7n;r77T. — , ,2 ■ (^2) 
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Appendix B: Fisher matrix for amplitude parameters 

Let us consider the gravitational-wave signal s of the form 



s(t;A) = ^Afe hk{t), 



(Bl) 



fc=i 



where the vector A collects the amplitude parameters, A := (Ai, A2, A^, and the known functions hk {k = 
1, . . . , 4) are given in Eqs. (jH). We further assume, as in Sec. 2, that the noise spectral density is constant (and equal 
to 5*0) over the bandwidth of the signal and that the approximations (jlOp are valid. Then the Fisher matrix for the 
signal's parameters A reads 



r(A) - ^ 



and its inverse is qual to 



r(A) 



5*0 



To{{a^){b^)-{aby 



f{a^) {ab) \ 

{ab) (62) 

(a2) (ab) 

V (ab) {b'^)J 



( {b^) -{ab) 
-{ab) (a2) 




(B2) 



V 







\ 



{ab) 
{ab) {a^)J 



(B3) 



Let us introduce new set of parameters 6 := {Hq, (j)Q,ip, i). Then the Fisher matrix T{6) for these parameters can be 
computed as (T denotes here matrix transposition) 



r(0) = j'^ ■ r(A) • J, 



(B4) 



where the Jacobi 4x4 matrix J has elements dAi/dOj — 1, ... ,4), which can be computed by means of Eqs. 
(El. 
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